Vertex cover problem studied by cavity method: Analytics and population dynamics 
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We study the vertex cover problem on finite connectivity random graphs by zero-temperature 
cavity method. The minimum vertex cover corresponds to the ground state(s) of a proposed Ising 
spin model. When the connectivity c > e = 2.718282, there is no state for this system as the 
reweighting parameter y, which takes a similar role as the inverse temperature /3 in conventional 
statistical physics, approaches infinity; consequently the ground state energy is obtained at a finite 
value of y when the free energy function attains its maximum value. The minimum vertex cover 
size at given c is estimated using population dynamics and compared with known rigorous bounds 
and numerical results. The backbone size is also calculated. 

PACS numbers; 75.10.Nr,89.75.-k,05.20.-y 
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I. INTRODUCTION 

The statistical physics of spin glass systems on 
infinitely connected lattices (e.g., the Sherrington- 
Kirkpatrick model) has been well understood and re- 
search interest is now focused on systems of finite con- 
nectivity (FC). For the limiting example of FC random 
lattice, development of the cavity method has been made 
in recent few years j3- In a FC random lattice, each 
vertex interacts only with a finite number of randomly 
chosen neighbours; and the local structure of the lattice 
is tree-like, with the shortest distance between two ran- 
domly chosen vertices diverging as the system size goes 
to infinity. This tree-like property makes it feasible to 
study the system by iterative cavity method The 
cavity method for n-dimensional regular lattice systems 
is yet to be worked out, and some aspects of such kind 
of systems has been understood using the well-developed 
mathematical tool of gauge transformation 0, |^ . 

The zero-temperature property of random lattice spin 
glasses is especially interesting. Here, the cavity formal- 
ism could be greatly simplified because only the mini- 
mum energy states are relevant Furthermore, many 
combinatorial optimisation problems in computer sci- 
ences could be studied through a mapping into an ap- 
propriate random spin glass model, examples of which 
include the K-sat JJ, the XOR-sat |3, , the vertex cov- 
ering problem 0,0] , the number partition problem |0 , 
and so on. Recently, the zero-temperature cavity method 
was used on the random K-sat problem and the phase di- 
agram for fc = 3 was obtained [Tsl Q . 

In this work we hope to improve our understanding on 
the vertex cover problem. Consider a random graph G 
composed of N vertices V = {1,2, - ■ ■ , TV}. Between any 
two vertices an edge is present with probability c/( — 1) 
and absent with probability 1 — c/{N — 1), so that on 
average each vertex has c neighbours (i.e., the average 
connectivity is c). The resulting edge set of graph G is 
denoted by E{G). A vertex cover of this graph consists 
of a set of vertices A = {ii, 12, • • • ,im}, with the property 
that if edge G E{G), then either i e A or j e A or 
both In the general case there are many different 



ways to cover a graph of size N; an interesting question 
is: Does there exist a vertex cover of size not exceeding 
xN {0<x< 1)7 

For large systems it was revealed that a sharp thresh- 
old value Xc{c) exists. When x > Xc{c), with probability 
approaching unity a vertex cover of size < xN could be 
constructed for a given graph; while when x < Xc(c) 
this probability approaches zero. This sharp threshold 
is also closely related to computational complexity. The 
vertex cover problem is NP-complete and a time 

growing exponentially with N may be needed to deter- 
mine whether a vertex cover of size < xN exists or not 
for a given graph In the case of FC random graphs, 
when X > Xc{c) or x < Xc{c), it is relatively easy for a 
heuristic algorithm to check the existence of such a ver- 
tex cover; however, when x ~ Xc{c), search complexity 
increases dramatically For the practical purpose of 
designing better algorithms, it is important for us to un- 
derstand the reason of this easy-hard-easy transition and 
to obtain a precise estimate of the threshold value Xc{c). 

An rigorous bound exists for xdc) 2;((c) < 

Xc{c) < 1 — In c/c, where xi(c) is the root of 



xi{c) lna;((c) 



- [1 - xi{c)]Hl - xi{c)] 
+ ic/2)[l ~ xi{c)]^ = 0; 



(1) 



furthermore, xJc) 
form at large c 



approaches the following asymptotic 



Xc{c) = 1 - (2/c)[lnc- lnlnc+ 1 - ln2] +o(l/c). (2) 

Using replica method of statistical physics, an analytical 
expression for Xc was found in p^ : 



Xcic) = 1 - W{c)/c - W'^{c)/2c, 



(3) 



where W{c) is the Lambert- W- function defined by 
W{c) exp[W{c)] = c. Equation @ is exact when c < 
e = 2.718282 ^J. For c > e Eq. © underestimates the 
true threshold value, and for c > 20.7 it is lower than the 
rigorous lower bound Eq. 

The present work focuses on the case of c > e. Using 
zero-temperature cavity method, we calculate both ana- 
lytically and numerically the value of Xc{c) and compare 
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it with Eqs. lU and (0) and with numerical calculations 
reported in p_Qi |. In Section^ an energy functional is in- 
troduced. In Section ITnl the cavity formalism is outlined 
and the free energy expression is given. We investigate 
the y — > oo situation in Section IIVI (y is a reweighting 
parameter jElHIIl, it plays the role of the inverse tem- 
perature [3 of conventional statistical mechanics). Section 
Ivlreports the population dynamics results on the thresh- 
old value Xc{c) and on the backbone size. We conclude 
the work in Section IVTI 



II. 



THE ENERGY FUNCTIONAL 



We attach to each vertex of the random graph G an 
Ising spin a — { — Associated with each spin 
micro-configuration is the following energy functional 



N 



A 



i=i (»j)es(G) 

where A is an constant parameter chosen to be greater 
than unity 0|. Denote N^-yc{G) as the size of the mini- 
mum vertex cover(s) of graph G, then the minimum en- 
ergy over all the 2^ possible spin configurations of graph 
G is 



c(G) - N. 



(5) 



Some explanation on Eq. First, E'min is reachable. 
We denote Ai„vc(G) as (one of) the minimum-sized ver- 
tex cover(s), and assign a = — 1 to vertices in Amvc(G) 
and cr = +1 to vertices outside. The energy of this spin 
configuration is 2Nyncv{G) — N. Second, no spin configu- 
ration could have lower energy. To see this, suppose that 
another spin configuration has lower energy than Eq. jSJ. 
This spin configuration must contain Na < N^^^ciG) neg- 
ative spins, with Na + XNq < iV,„vc(G), where A^o = 
(1/4) Z](i.j)6_E(G)(l + '^0(1 + o-j)- However, one at most 
need to change the spin values of Nq vertices from -f 1 to 
— 1 to make the sum j)£E{G)(^^'^'')i^^'^j) ^ ^^'^ 
the resulting new set of negative spin vertices is a vertex 
cover with size Na + Nq < Ai„vc(G). This conflicts with 
our original assumption that Anwc(G) is the minimum 
vertex cover size. 

The problem of finding Xc{c) is converted to finding 
the average of -Emin/-^ over the random graphs G; 



Xaic) = il/2){l+E^in/N). 



(6) 



Here (•) means the average of (•) over different realiza- 
tions of the random graph. We use cavity method to 
estimate E^i^/N. In the next section, the cavity formal- 
ism for the present problem is outlined. The reader is 
referred to 0, 0, 0, ITsL Q| for more detailed discussion. 



III. THE ZERO-TEMPERATURE CAVITY 
FORMALISM 

At zero temperature, only the minimum energy config- 
urations are relevant. There could be a great many en- 
ergy local-minima for Eq. Q . For very large system size 
N, we group these configurations into different "states" . 
A state of the system corresponds to a set of spin micro- 
configurations. These spin configurations all have the 
same energy, which is a local minimum of Eq. |0J; and 
two such spin configurations differ only in a finite num- 
ber of spin flips. The average number of states at given 
density e of local minimum energy is assumed to be an 
exponentially increasing function of system size A'", and 
is characterised by the entropy density I](c, e). We can 
introduce a reweighting parameter y and define an zero- 
temperature free energy density ^{y) through the follow- 
ing equation 



deexp[-Nye + NY.{c,e)] = exp[-A^y$(y)]. (7) 



Equation ||7J) has the same form as the conventional defi- 
nition of free energy in textbooks of equilibrium statisti- 
cal physics. The reweighting parameter y plays the role 
of inverse temperature. A large value of y ensures that 
states with lower energies will be favoured, provided that 
such states exist (i.e., S(c, e) > 0). 

Suppose we have a system of N vertices (spins). Now 
add a new spin cto into the system and connect it to k 
preexisting spins cti , • • • , (Jfc , where k obeys the Poisson 
distribution of mean c, Pp{k,c) — e'~'^c^ /k\. The energy 
of the A^-spin system at fixed value of the spins cti , • • • , cTfc 
is supposed to be 



i?W( 



CTi, 



,fTfe) =^ 



k 

E 



(8) 



with A being a constant. In the above equation, hi is 
the local field (called the cavity field) felt by spin ai in 
the absence of do in a given macroscopic state. Since the 
graph is locally tree-like, as the system size becomes very 
large, the shortest distance between two randomly chosen 
cavity spins also becomes large; therefore, the cavity field 
hi felt by spin ai becomes independent of the values of 
all the other cavity spins . After the addition of spin 
(To, the minimum energy of the (A -I- l)-spin system at 
fixed value of is 



A - (To, 



if = 



Here, 



w{h) =0 if /i = 1 and = \h\ if < 0, 
u{h) = -1 if /i = 1 and = li h < 0. 



(9) 



(10) 



(We have used the fact that the cavity fields at zero tem- 
perature are integer- valued and do not exceed unity |21|.) 



3 



The energy shift caused by the addition of spin (Tq is 
AEi = -1 (if fc = 0) and AEi = -Y^Lii^if^i) ~ ~ 
\^ + J:■=lHh^)\ (iffc>l). 
Equation (jnj indicates that the cavity field at spin cto 

is /lo = 1 (if fc = 0) or ho = 1 + (if ^ ^ !)■ % 

definition |21| , the cavity field of a spin at each state has 
unique value; however, its value may be different for dif- 
ferent states. Denote the probability distribution of the 



cavity field at vertex i amon g dif ferent states as Pi{h) (it 
is called the h-survey in Ref. [3). Because different ver- 
tices have different local environments, the h-surveys are 
different for different vertices. With the introduction of 
the reweighting parameter y which favours lower-energy 
states, the h-survey at spin ctq is related to those of the 
cavity spins by 



Po(/j) = '5°<5(/^ - 1) + [1 - <5^]C / ll[E,ih,)dh,]S[h - 1 -Y,Hh^)]eM-y^E,), (11) 

i=l 1=1 

where C is an normalisation constant. Eq. Hll() is actually a self-consistent equation for the h-survey. A careful 
analysis of Eq. (|ll|l leads to the following expression 

C/(5(/i + with probability pi 

P{h) ^ { 5{h- 1), with probability p2 (12) 

a X];^Q C/(5(/i + Z) + (1 — a)5{h — 1), with probability py, = 1 — pi — P2 



r 



where > and J2i C/ = 1; " £ (0, 1) is determined by 
certain probability distribution p{a), and 

pi = 1- exp(-c]32), 

P2 = exp(-c(l -pi)). (13) 

For very large system size N, Eq. suggests that 



at fixed value of y, $(y) = e — E(c, e)/y, with e being 
implicitly determined by 9I](c, e)/9e — y. An explicit 
expression for the free energy density could be obtained 
by the following way: 

After the addition of spin ctq, the averaged number of 
states of the {N + l)-system is 



where P'^^^(Ai?i) is the probability distribution function of the energy shift AEi: 

p(i)(Ai?i) = <5^<5(Ai?i + 1) + (1 - 5'1) J Y[[n{h,)dh,]6[AEi + ^(u;(/i,) - \h,\) + |1 + ^u(/i,)|]. (15) 

i—l i—1 i—1 

A logarithm operation is performed on Eq. (|14|) , and the resulting equation is averaged over all the possible realizations 
of the cavity fields and k (this average operation is denoted by an overbar in the following equation). We obtain that 



-2/$(2/) = E(c, e) - 2/e = -c^^^ + In J dAEiP(^) {AEi) cxp{-yAEi). (16) 

To compute dT,/dc, we setup an edge between two cavity spins. The energy shift, AE2, caused by this new edge 
obeys the following distribution 

r 2 ^ 
P^^\AE2)^ / [][d/i,P,(/i,0]'5[AS2-(min,,,,,[-(l + fTi)(l + (72)-/^iai-/i2a2] + |/ii| + |/i2|)]. (17) 

i=i 

The averaged number of states of the new system is 
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After performing the same procedure as mentioned below Eq. (|15|) . we find that 

9I](c,e) 



dc 



(l/2)hi J dAE2P^^KAE2)eM-y^E2). 



(19) 



The free energy density expression could be obtained from Eqs. H16|) and (|19|l . After taking into consideration Eq. H12|l , 
we arrive at the following expression for the free energy density: 



$(y) = 2pi - 1 - cpl + / dap{a) ln(a + (1 - a)e-^y) 



P3 

y 



E^^/nW«)<i«lln(e-' + (l 



m—l 
2 - 2 



i=l 



/ flida^pia.)] HI - (1 - e-'y) f[il - a,)). 



(20) 



At given y, the energy density and the entropy density are calculated by e = <^{y) ^ yd<^{y) / dy and E — y^d^{y)/dy, 
respectively. 



IV. THE y^oo LIMIT 



or 



At this stage, it is helpful for us to introduce an auxil- 
iary probability distribution function called the u-survey 



g,(M) = C / dhP{h)S{u~u{h))exp{y{w{h)-\h\)). (21) 



g,(w) = f]S{u) + [1 - r]]S{u + 1). 



(22) 



The hybrid h-survey and u-survey at a given vertex are 
Based on Eq. H12|l we know that Qi{u) = 5{u) (with related by a = 77/(77 -I- (1 — 77)6^). The distribution of the 
probability pi) or Qi(u) = 5{u+\) (with probability P2), V value in the hybrid u-survey Eq. H22|l is governed by 

I 



m=l •' i=l ri 



1=1 



Uh + ey{i~vi)] + iey ~i)Um 



-)■ 



(23) 



The free energy expression of Eq. 1)20(1 could be rewritten in the following form 



rf.^ \ ^2 , / n2i , P3(1+CP2) 

*(j/) = Pi -P2 - 7:[(1 -Pi) +{P2) ]H 



ln[7; 4- e-y{l - 77)] 



P3[l + c(l-pi)] 



2 

ln[l - 77 + 6-^77] + ^ln[e-y + (1 - e-!')(77i + 7/2 - 27717^2)]. 

2y 



r 



(24) 



The overbars in Eq. H24() denote the average over the 77 
distribution given by Eq. ((23|l . 

When c < e, the only solution of Eq. 1(1 3|l is that pi — 
1 — P2, P2 — exp(— CP2) — W{c)/c, and ps ~ 0. In this 
case we recover Eq. Q. S(c, e) — for this solution, 
indicating that there is only one state. It was found that 
for c < e, all the vertices of the random graph could be 
removed by application of an leaf- removal algorithm |l9j | . 



When c > e the above-mentioned solution becomes 
unstable, and a new solution appears with p2 — 
exp[-cexp(-cp2)] < exp(-cp2), Pi = 1 - exp(-ciP^) 
and p3 > 0. The analysis of Bauer and Golinelli [13 
reveals that a core with size proportional to N remains 
after application of the leaf-removal algorithm. The core 
is a strongly connected subgraph, in which each vertex is 
connected at least to two other vertices. 



FIG. 1: The r] value distribution Eq. I|28|l obtained by pop- 
ulation dynamics at c = 6.0. The p(r;) distribution strongly 
depends on the value of the reweighting parameter y. In this 
figure, three curves at y = O.f (dotted line), y = y* = 1.5545 
(corresponding to the free energy maximum, dashed line), and 
y — 20 (solid line) are shown. 



To estimate the value of Xc{c) for c > e, let we first 
consider the Hmiting situation of j/ ^ oo. Equation {T)) 
ensures that the minimum e corresponds to y = cx), pro- 
vided that the configurational entropy is non-negative at 
this Umit. 

From Eq. (|23() we know that as j/ ^ cxo, 

p{7]) = rid{7] - 0+) + r2S{7] - 1-) + r3p*{i]), (25) 

where p*(r]) is a uniform distribution over (0, 1). Equa- 
tion (|25|l is confirmed by the population dynamics calcu- 
lation (Fig.^l. It is easy to verify that 

1 — cp2 c + cpi — In c In c — 1 

n = , r2 = , rs = . 

CP3 cp3 cpz 

(26) 

Consequently, we obtain from Eq. 124|) that 

too = 1 - hi^ (c) /2c - ln(c) /c - 3 /2c, 
Soo = -7r2(lnc-l)Vl6c. (27) 

The minimum energy value given in the above equation 
is an improved lower-bound of the true average minimum 
energy. It was obtained also in by replica method. 
However, at y = cxd the entropy density Sqo is negative, 
suggesting that there is no state at this energy density. 
Indeed the Xc[c) value [Eq. lO] of this solution also ex- 
ceeds the rigorous lower-bound when c > 27.3. The true 
energy density must be higher than the value given in 
Eq. H27I) . A better estimate of the minimum energy den- 
sity could be obtained by calculating the maximum value 
of $(y) with respect to the reweighting parameter y j2^ . 
This is done in the next section with population dynam- 
ics 0. 



FIG. 2: The minimum vertex cover size Xc(c) estimated by 
population dynamics (solid line and squares) and its com- 
parison with the exact numerical enumeration values of 
(circles), the replica-symmetric estimation Eq. (|^ (dotted 
line), and the estimation given by Eq. H27|l (dashed line). In 
the inset, the Xc{c) value obtained by population dynamics 
(solid line) is compared with the rigorous lower-bound Eq. Q 
(dashed line) and the asymptotic value Eq. Q (long dashed 
line) . 



TABLE I: The Xc(c) value obtained by the population dy- 
namics (column 2) and its comparison with the exact nu- 
merical enumeration value reported in [l^ (column 4). The 
reweighting parameter value y* of column 3 at given c corre- 
sponds to the maximum of the free energy density Eq. 1)24^ . 
To determine y*, Eq. 124^ is fitted with two adjustable pa- 
rameters: ^(j/) = too + C\ly + C2 exp(— ?/)/i/, where too is 
given in Eq. H27|l . The population size adopted in the present 
work is 20, 000. 



c 


Xc (cavity) 




(enumeration) 


4.0 


.5194 


1.3093 


.523 ± .001 


5.0 


.5603 


1.5758 




6.0 


.5934 


1.5545 


.599 ± .001 


7.0 


.6210 


1.5260 




8.0 


.6443 


1.5259 


.656 ± .003 


9.0 


.6643 


1.5326 




10.0 


.6819 


1.5434 


.697 ± .003 



V. POPULATION DYNAMICS AT FINITE y 

To obtain the value of $(y) at any given y, the tech- 
nique of population dynamics is used 2] . A large popula- 
tion of rf'a is generated and this population then evolves 
according to Eq. (|23|) . The averages in Eq. (|24|l are cal- 
culated numerically. The resulting estimates of Xc(c) are 
shown in Fig. |21 and listed in Table |l| 

The threshold curve Xc(c) obtained by the present 
method lies within the rigorous bound given by Eq. 
It is therefore an improved estimate compared with 
Eq. Q and Eq. (|27|) . both of which exceed the rigorous 
lower-bound when c is larger than certain value. How- 
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FIG. 3: Fraction of frozen spins calculated by population 
dynamics j23| (squares) and its comparison with numerical 
result reported in pl| and the replica-symmetric result 
(thick solid line). 

ever, the values of Xc estimated by the cavity method 
is systematically smaller than enumeration results on fi- 
nite systems (Fig.|51). When the average connectivity c is 
large, the estimated Xc(c) value approaches the asymp- 
totic value but it lies below the asymptotic curve (Fig. [21 
inset). These discrepancies suggest that the Xc(c) thresh- 
old value obtained by the cavity method is not exact; 
it could serve as an improved lower-bound of the real 
threshold curve. The reason for the failure of the cavity 
method to obtain exact threshold values for c > e are 
discussed in the next section. 

The cavity field distribution Eq. H12|l contains infor- 
mation on the spin state at a given vertex. A spin 
CTi will be fixed at ai ~ +1 if the local cavity field is 
distributed according to P{hi) = S{h — 1) or P{hi) = 
aEi=o O^i^ + + (1 - - 1) but with a ^ 0+ 

[23]; on the other hand, will be fixed at = —1 
with probability J^'^iCi if P{hi) = Sz^oO^Ci + or 
P{hi) = ay2'ZoCi5{h + 1) + {\ - a)S{h - 1) but with 
a 1^ 23]. For the minimum vertex covers at y = y*, 
the probability for a randomly chosen vertex to have 
fixed spin value is calculated by population dynamics (see 
Fig. |3J). The fraction of frozen spins calculated by the 
population dynamics is much lower that that obtained 
by numerical enumeration This discrepancy may 

be further indication that the full hierarchy of replica 
symmetry breaking is needed to completely describe the 
property of the minimum vertex covers. We have noticed 
that the fraction of frozen spins strongly depends on the 



value of the reweighting parameter y. When y ^ oo the 
replica-symmetric result of [TTI | is recovered; while for 
y ~ 2y*Sy* the data obtained by numerical enumera- 
tion is approached. 

VI. DISCUSSION 

In this work, we estimate the average minimum vertex 
cover size Xc{c) for random graphs of finite connectivity 
c > e in the large N limit. The obtained Xc{c) curve 
lies within the rigorous bound |l7j and approaches the 
asymptotic curve Eq. (|2Jl at large c values. It could be 
regarded as an improved lower-bound of the real thresh- 
old Xc{c) value. 

When c > e the threshold Xc estimated by the cav- 
ity method is not exact. The reason may be the follow- 
ing: The cavity method is equivalent to one-step replica- 
symmetry-breaking iSi »SJ 7 ^-nd the cavity fields on dif- 
ferent vertices are considered as uncorrelated. Because 
of the core percolation beyond c = e Toj in the random 
graph, the cavity fields of different vertices may actually 
be correlated strongly. To partly account for this effect, 
one possibility is to consider also non-integer cavity fields. 
We hope to return to this point in a latter work. 

The cavity method has inspired very efhcient algo- 
rithms to tackle the random K-sat problem. In the 
random K-sat problem, there exists a glassy phase at 
y ^ oo for certain range of the connectivity c 0, . 
In this glassy phase, the minimum energy is still located 
at y — oo but the complexity is positive, and the algo- 
rithm based on the idea of cavity field 1a\ works well 
in this phase. For the vertex cover problem, the present 
work suggests that such a y — oo glass phase does not 
exist (this conclusion seems also true for the vertex cover 
problem on random hyper-graphs where each "edge" is 
a triangle This indicates that vertex covers of 

size Nxc{c) + 0{1) are extremely few. It remains to be 
seen whether or not algorithms based on cavity method 
could efficiently find vertex covers of size slightly beyond 
Nxc{c) for a given random graph. 
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